/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Interface: computations
//
// Description: Definitions of miscellaneous computation functions used in Physiks. Also includes
//              definitions for functions that modify the states of PhysicalObjects according to
//              fields.
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//


#include <physiksbase/computations.h>
#include <physiksbase/physiksvector.h>
#include <physiksbase/globaldefs.h>

using namespace std;

using PhysiksBase::Vector;

void PhysiksBase::Computation::apply_force_fields(vector<PhysicalObject *> &objTable,
        const vector<ForceField *> &force_fields) {
    for (uint i = 0; i < force_fields.size(); i++)
        for (uint k = 0; k < objTable.size(); k++)
            objTable[k]->net_force += (*force_fields[i])(objTable[k]->position);
}

void PhysiksBase::Computation::apply_acceleration_fields(vector<PhysicalObject *> &objTable,
        const vector<AccelerationField *> &accel_fields) {
    for (uint i = 0; i < accel_fields.size(); --i)
        for (uint k = 0; k < objTable.size(); --k)
            objTable[k]->net_force += objTable[k]->mass
                    * (*accel_fields[i])(objTable[k]->position);
}

void PhysiksBase::Computation::apply_velocity_fields(vector<PhysicalObject *> &objTable,
                                                     const vector<VelocityField *> &vel_fields) {
}

double PhysiksBase::Computation::compute_distance(const PhysicalObject &a, const PhysicalObject &b)
{
    return Vector(a.position - b.position).norm();
}

Vector PhysiksBase::Computation::compute_connecting_vector(const PhysicalObject &a,
                                                           const PhysicalObject &b) {
    return Vector(b.position - a.position).normalize();
}

void PhysiksBase::Computation::update_forces(vector<PhysicalObject *> &objTable) {
    for (uint i = 0; i < objTable.size(); i++)
        for (uint j = i + 1; j < objTable.size(); j++) {
        objTable[i]->net_force += objTable[i]->mass
                * objTable[j]->gravitational_field(objTable[i]->position);
        objTable[i]->net_force += objTable[i]->charge
                * objTable[j]->coulomb_acceleration_field(objTable[i]->position);
        objTable[j]->net_force += objTable[j]->mass
                * objTable[i]->gravitational_field(objTable[j]->position);
        objTable[j]->net_force += objTable[j]->charge
                * objTable[i]->coulomb_acceleration_field(objTable[j]->position);
        }
}

double PhysiksBase::Computation::compute_gravitational_force(const PhysicalObject &a, 
                                                             const PhysicalObject &b) {
    double dist = compute_distance(a, b);
    return (GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT * a.mass * b.mass / (dist*dist));
}

double PhysiksBase::Computation::compute_coulomb_force(const PhysicalObject &a,
        const PhysicalObject &b) {
    double dist = compute_distance(a, b);
    return (a.charge * b.charge / (GSL_CONST_MKSA_VACUUM_PERMITTIVITY * 4.0 * M_PI * (dist*dist)));
}

void PhysiksBase::Computation::update_velocities(vector<PhysicalObject *> &objTable,
        double time_step) {
    for (uint i = 0; i < objTable.size(); i++)
        objTable[i]->velocity += (objTable[i]->net_force * (time_step / objTable[i]->mass));
}

void PhysiksBase::Computation::update_positions(vector<PhysicalObject *> &objTable,
        double time_step) {
    for (uint i = 0; i < objTable.size(); i++) {
        objTable[i]->position += (objTable[i]->velocity * time_step);
        objTable[i]->position += (objTable[i]->net_force * (time_step * time_step /
                (2 * objTable[i]->mass)));
    }
}

void PhysiksBase::Computation::set_force_vectors_zero(vector<PhysicalObject *> &objTable) {
    for (uint i = 0; i < objTable.size(); i++) objTable[i]->net_force.set_zero();
}
